Search-and-rescue in the Central Mediterranean Route does not induce migration: Predictive modeling to answer causal queries in migration research

State- and private-led search-and-rescue are hypothesized to foster irregular migration (and thereby migrant fatalities) by altering the decision calculus associated with the journey. We here investigate this ‘pull factor’ claim by focusing on the Central Mediterranean route, the most frequented and deadly irregular migration route towards Europe during the past decade. Based on three intervention periods—(1) state-led Mare Nostrum, (2) private-led search-and-rescue, and (3) coordinated pushbacks by the Libyan Coast Guard—which correspond to substantial changes in laws, policies, and practices of search-and-rescue in the Mediterranean, we are able to test the ‘pull factor’ claim by employing an innovative machine learning method in combination with causal inference. We employ a Bayesian structural time-series model to estimate the effects of these three intervention periods on the migration flow as measured by crossing attempts (i.e., time-series aggregate counts of arrivals, pushbacks, and deaths), adjusting for various known drivers of irregular migration. We combine multiple sources of traditional and non-traditional data to build a synthetic, predicted counterfactual flow. Results show that our predictive modeling approach accurately captures the behavior of the target time-series during the various pre-intervention periods of interest. A comparison of the observed and predicted counterfactual time-series in the post-intervention periods suggest that pushback policies did affect the migration flow, but that the search-and-rescue periods did not yield a discernible difference between the observed and the predicted counterfactual number of crossing attempts. Hence we do not find support for search-and-rescue as a driver of irregular migration. In general, this modeling approach lends itself to forecasting migration flows with the goal of answering causal queries in migration research.


Data
presents the variables used to construct our target time series of attempted crossings, as well as the covariates used in our prediction model, together with their sources and operationalization.

Dates of search-and-rescue NGO-led operations and EU-led operations
Our information on the dates of search-and-rescue NGO missions taking place -including the days when a rescue mission started and when vessels were being held in a harbor or ending a mission -is collected from multiple sources. The first source of information was the official websites of the NGOs, where sometimes this information was made available. However, given that only few of them provided this information, we made used of published news prints and research articles containing the exact days. This helped to get the schedules of the different search-and-rescue missions. Although at times these strategies did not provide exact dates, they gave us with at least a rough estimate of the time frame. This estimate was further used as a reference period to further search in other sources. Many of the NGOs engaged in search-and-rescue missions make use of social media platforms, like Twitter, to raise awareness of the events occurring in the Mediterranean ocean, inform their supporters of their activities, or just collect funds for future rescue attempts. Given that based on the estimated time frame, we knew they were deployed in the central Mediterranean Sea, we made use of the information in social media platforms as a second direct source of information. If it was not possible to find dates through any official sources the research had to be widened to include secondary more indirect sources like interviews with NGO officials or reports in newspapers. Together these three different sources provided the exact dates of start, end and (involuntary) breaks in the several search-and-rescue missions in the Mediterranean Ocean. Binary indicators taking the value of one when the NGO was active in the Central Mediterranean and zero otherwise.

Deads and missing count
One important factor affecting our measure of attempted crossings is the underregister of deaths. The true number of deaths in the CMR is an unknown quantity that may affect our estimates 1 . In Figure S3, we show the evolution of the ratio of the sum of known deaths and missing migrants over the sum of arrivals and pushbacks along the CMR for the period of study. This figure suggest that deaths are a relatively small component of our measure of attempted crossings, particularly so during our intervention periods. In a recent study 2 , underregister was found to be of a relatively small magnitude: less than 10% of deaths are estimate to have gone unrecorded. It is, however, unclear how this underregister changes over time, which is crucial to evaluate how it might affect our estimates.
It might well be that underregister is higher during the intervention periods under study, but the data in Figure S6 suggest that before search-and-rescue activities the ratio of deaths and missing was far higher than during search-and-rescue activities, when the ratio seems to be in fact much lower and around 0.1. Equally so after the start of the cooperation between the Libyan Coast Guard and the European Union. There is reason to believe that search-and-rescue activities in fact could reduce underregister of deaths given that more efforts are made to search for fatalities, thus improving reporting in a more systematic fashion. Hence, mostly the pre-intervention period would be affected by underregister. We believe this is unlikely to affect our estimates given that our BSTS models would be able to correctly capture the pre-intervention dynamics with the true unobserved number of deaths.
Similarly to the main study on the effects of the three intervention periods -The Mare Nostrum period, the private-led search-and-rescue, and the period of

Model validation and comparison with other machine learning models
Our causal identification strategy is based on a combination of the difference-in-differences and synthetic counterfactual frameworks 4,5 , which, as argued in the main paper, is the most suited model for the type of question we try to answer given the data and intervention settings. In our case, where we have pre and post intervention periods with different regimes of search-and-rescue and pushbacks being active/inactive, we want to be able to establish a fair comparison while attending to basic elements of aggregate, discrete count time series data used for this (i.e., cycles or changes in variability, trends, seasonality, and a randomness).
The results of our BSTS models depend to an important extent on which state specification is used 5 , and also, though to a lesser extent, on the selection of priors for many parameters. Model validation is therefore important for these state specifications. In particular, we perform model validation to evaluate the performance of these main components of the state specification in our BSTS models: • The local level given by the equation α t+1 = α t + ϵ t , with ϵ t ∼ N (0, σ); and where cross-validation is applied on a set of parameters for σ.
• The local linear trend component which assumes mean and slope of the trend follow a random walk, with equations for the mean given by µ t+1 = µ t + δ t + ϵ t , with ϵ t ∼ N (0, σ µ ), and for the slope δ t+1 = δ t + η t , with η t ∼ N (0, σ δ ).
• The semi-local linear trend, though equally assuming the same random walk process for the mean, assumes for the slope a stationary AR1 process centered on a value D, as expressed in , which is particularly useful when making projections far into the future.
In hindsight, adding the seasonality component to the model made the fitting worse given that the seasonality, although present in the data, is not stationary, meaning that it changes over the period, and which correspond to an important feature of the target series to be predicted by our models. Therefore, we excluded this component from the final model.
We perform cross-validation using 5 folds and evaluate the prediction employing a different horizon length depending on the availability of pre-intervention data (three months for the Mare Nostrum and six months for the private-led search-and-rescue and the extension of the Libyan search-and-rescue area). Given that we have three intervention periods and hence three different BSTS models, we perform model validation for each of these models employing only the pre-intervention data in each case. We evaluate the performance of the different parameter configurations employing three metrics: root mean squared error (RSME), the mean absolute prediction error (MAPE), and the mean absolute scaled error (MASE).  The comparison of cumulative absolute error shown in Figure S6 reveals that the default model beats the cross-validated measures for the components of the state space model we focused on. This behavior is probably due to the short horizon on which we can evaluate the performance of the model. Our predictions for a causal effect extend over years and therefore the fit of the model should be evaluated on a much longer horizon. This, however, is not possible due to the limited availability of data. Estimates of the different models, however, do not differ much as shown in Figure S7. The most notorious difference between these estimates are the standard errors. Some models provide much wider CIs, but the effect is of a similar magnitude across models. The final models we present in the main paper have a relatively good fit as seen in the analysis of residuals presented in Figures S8 and S9. The autocorrelation plot in Figure S8 shows the characteristic drop towards zero.
Although the tails of the error distribution do not follow precisely on the straight line -probably due to Figure S8: Autocovariance/autocorrelation function for the predicted series for the three intervention periods variable omission, our models show a good fit in general and not a strong deviation from the assumed normal distribution for the residuals.

Selected covariates in BSTS final models
By means of spike-and-slab prior, for our three predictive models we have a subset of covariates for each intervention period which were deemed predictive. These are shown in Figure S3. The color of the bars representing the direction of the effect. For illustration, we present covariates with an inclusion probability of 10% or more, but other covariates with a lower inclusion probability and hence a lower weight were part of the model.

Mediation analysis
One of the most important assumptions of our approach is that the control series used in the BSTS model are exogenous with respect to the interventions of interest. If these control series were affected by any of our interventions, our identification strategy would be compromised, as we would be overcontrolling for potentially relevant elements of the causal process connecting changes in search-and-rescue politics and number of attempted crossings. In this section, we pay close attention to this issue by highlighting the potential for the labor market indicators (e.g., unemployment rates in Europe and job searches in North African countries) to be affected by the hypothetical "pull-factor" effect of search and rescue operations or the extension of the Libyan search-and-rescue area. We do this only for these variables because we believe these are the only ones we can plausibly see as potentially affected by the interventions in a feedback way (i.e., if search-and-rescue does constitute a "pull factor," more migrants will leave their country of origin or destination, reducing the labor supply in those labor markets, and more migrants will arrive to Europe, increasing the labor supply in corresponding local European labor markets). However, we do not believe the number of arrivals is large enough to have such disruptive effects on European or North African labor markets. In addition, there are physical and bureaucratic labor market access constraints for migrants entering Europe and seeking asylum.
In order to rule out this hypothetical mechanism, we exclude the labor market indicators from the set of control series out of which the BSTS model can select and re-run our models. The results of these models, shown in Figure S11, show very similar effects, as those including these series with only small variations in the predicted counterfactual, which makes intuitive sense given that each variable in our model has a minor contribution to the generation of the synthetic counterfactual.
Second, and to rule out potential mediators, we perform a series of bivariate Granger causality tests, which  Tables S4-S6 for the three intervention periods. For each period, each model contained hundreds of variables, but only for a few of these control series we were able to detect statistically significant associations with the rejection of the reverse causal hypothesis. Given that it is not plausible that the interventions have affected the behavior of these control series, we interpret these few significant values as spurious correlations.       Authors find no clear time trend in the overall mortality rate, and they find no evidence of a higher number of arrivals in the low search-and-rescue period than in the high period, as is expected from the "pull factor" claim.  The cut back on state-led search and rescue led to an increase in the mortality risk of crossing the CMR.